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Abstract 

We provide a priori error estimates for the spectral and pseudospectral 
Fourier (also called planewave) discretizations of the periodic Thomas- 
Fermi-von Weizsacker (TFW) model and for the spectral discretization 
of the Kohn-Sham model, within the local density approximation (LDA). 
These models allow to compute approximations of the ground state energy 
and density of molecular systems in the condensed phase. The TFW 
model is stricly convex with respect to the electronic density, and allows 
for a comprehensive analysis. This is not the case for the Kohn-Sham LDA 
model, for which the uniqueness of the ground state electronic density 
is not guaranteed. Under a coercivity assumption on the second order 
optimality condition, we prove that for large enough energy cut-ofFs, the 
discrctizod Kohn-Sham LDA problem has a minimizer in the vicinity of 
any Kohn-Sham ground state, and that this minimizer is unique up to 
unitary transform. We then derive optimal a priori error estimates for 
the spectral discretization method. 



1 Introduction 

Density Functional Theory (DFT) is a powerful method for computing ground 
state electronic energies and densities in quantum chemistry, materials science, 
molecular biology and nanosciences. The models originating from DFT can 

be classified into two categories: the orbital-free models and the Kohn-Sham 
models. The Thomas-Fermi-von Weizsacker (TFW) model falls into the first 
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category. It is not very much used in practice, but is interesting from a math- 
ematical viewpoint [1, 7, 12]. It indeed serves as a toy model for the analysis 
of the more complex electronic structure models routinely used by Physicists 
and Chemists. At the other extremity of the spectrum, the Kohn-Sham models 
[8, 11] are among the most widely used models in Physics and Chemistry, but 
are much more difficult to deal with. We focus here on the numerical analysis 
of the TFW model on the one hand, and of the Kohn-Sham model, within the 
local density approximation (LDA), on the other hand. More precisely, we are 
interested in the spectral and pseudospectral Fourier, more commonly called 
planewavc, discretizations of the periodic versions of these two models. In this 
context, the simulation domain, sometimes referred to as the supercell, is the 
unit cell of some periodic lattice of M^. In the TFW model, periodic boundary 
conditions (PBC) arc imposed to the density; in the Kohn-Sham framework, 
they are imposed to the Kohn-Sham orbitals (Born- von Karman PBC). Impos- 
ing PBC at the boundary of the simulation cell is a standard method to compute 
condensed phase properties with a limited number of atoms in the simulation 
cell, hence at a moderate computational cost. 

This article is organized as follows. In Section 2, we briefly introduce the func- 
tional setting used in the formulation and the analysis of the planewave dis- 
cretization of orbital-free and Kohn-Sham models. In Section 3, we provide a 
priori error estimates for the planewave discretization of the TFW model, in- 
cluding numerical integration. In Section 4, we deal with the Kohn-Sham LDA 
model. 



2 Basic Fourier analysis for planewave discretiza- 
tion methods 

Throughout this article, we denote by T the simulation cell, by TZ the periodic 
lattice, and by TZ* the dual lattice. For simplicity, we assume that T = [0, L)^ 
{L > 0), in which case U is the cubic lattice LZ^, and Tl* = "^1? . Our 
arguments can be easily extended to the general case. For fc € 7?.*, we denote 
by efe(x) = |r|~^/^ e^^'^ the planewave with wavevector k. The family {ek)keTZ' 
forms an orthonormal basis of 



L|(r,c) 



{u G L^oc(M^C) I u 7^-periodic} , 



and for all u G L'i 



(r,c). 





u{x)e 



— ik-x 



dx. 
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In our analysis, wc will mainly consider real valued functions. We therefore 
introduce the Sobolev spaces of real valued 7?.-periodic functions 

H^{T) := luix) = J2^kekix)\ (1 + < oo and Vfc, = n*^ I , 

I k£TZ* k£n* ) 

s e M (here and in the sequel a* denotes the complex conjugate of the complex 
number a), endowed with the inner products 

For A^c € N, we denote by 



[ken-\\k\<^N, ) 

(the constraints C-k = c% imply that the functions of Vn^ are real valued). For 
all s e K, and each v G H^{T), the best approximation of v in Vn^ for any 
H^-norm, r < s, is 

ken" I |fc|<^JVc 

The more regular v (the regularity being measured in terms of the Sobolev 
norms H'^), the faster the convergence of this truncated series to v: for all real 
numbers r and s with r < s, we have for each v € H^(r), 

^ {ij '^c-^^-^^ii-iihj. (2) 

For Ng g N \ {0}, we denote by ^^'^'-'^3 the discrete Fourier transform on the 
carterisan grid GNg '■= Wg'^ °^ function ^ e C^(r,C), where 

C^(r,C) := {u e C°(R^C) I u 7e-periodic} . 

Recall that if <^ = "Yliken-" ^k € C^(r, C), the discrete Fourier transform of <p 
is the A'g7?.*-periodic sequence ^^'^'^9 = {(t>^^'^'^'')ke'R' where 



7-FFT,JV„ 

=jp 2^ nx)e -- = |i| -'- 2^ (Pk+NgK. 



We now introduce the subspaces 



Spanje^'^ I / G ^Z, i/| < ^ (\^)] (^9 odd), 

Span j e^'^ M € ^Z, |;| < ^ ( j © C(e*-^»^/^ + e-^^^««/^) {Ng even), 
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e C^([0, L), C) and d\m{W}P^) = N,), and WfP^ = WjJ^ ® W}^^ ® W}P^. 
Note that is a subspace of i?|t(r,C) of dimension A^^, for all s e M, and 
that if Ng is odd, 

= Span |e, I fc e 7e* = ^Z^, |fcU < ^ (^V^) } o^d). 

It is then possible to define the interpolation projector IjVg from C^(r,C) onto 
W|r° by [Jjv,(?i)](a;) = (i){x) for all x e Gn,- It holds 

The coefficients of the expansion of Xisi^{(p) in the canonical basis of WpP is 
given by the discrete Fourier transform of (p. In particular, when Ng is odd, we 
have the simple relation 

lMM) = n^' E ^r^'^^'ek {Ng odd). 

feeTC*||fe|oo<¥(^) 
It is easy to check that if (p is real- valued, then so is lNg{4>)- 
We will assume in the sequel that Ng > 4Nc + 1. We will then have for all 

[ V4N,= E i^l ^4JVe(a;)= / Xjv,(u4JvJ. (4) 

The following lemma gathers some technical results which will be useful for 
the numerical analysis of the planewave discretization of orbital-free and Kohn- 
Sham models. 

Lemma 2.1 Let N^ e N* and Ng € W such that Ng > AN^ + 1. 

1. Let V he a function o/C^(r,C) andvN^ andw^c he two functions ofVN^. 
Then 



j iNgiVVN^WNj = j lNg{V)VN^WN^\ (5) 

TNg{V\vNf) 

r 



< \\V\\lAvnA\1.. (6) 

# 



2. Let s > 3/2, < r < s, and V a function of H^{T). Then, 

\\T12nAInAV))\\^,^ < I^J^In^{\V\^)^ ' ■ (8) 

||n2iv.(x,v,(v^))L. < (i + c«,,)||y||H^, (9) 
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for constants Cr,s independent ofV. Besides if there exists m > 3 and C e 
M+ such that \ Vk\ < C\k\~"^, then there exists a constant Cy independent 



of Nc and Ng such that 



(10) 



3. Let (j) be a Borel function from ^+ to M such that there exists G M+ 
for which < C^{1 + 1^) for all t G M+. Then, for all vn, G Vn,_, 

I^InMM)) < (|r| + ll^ivJI^J . (11) 



Proof For Z2n^ € , it holds 



/ In,{Vz2nJ = E 

E 



L 



V{x)Z2N^{x) 



{InAV)){x)z2nAx) 



-L 



^N„ (V) Z2Nc 



(12) 



since IjVg(^)-22JVc € V^g+'iN^ C V2JVg is exactly integrated. The function wat^ujjv^ 
being in V2n^-, (5) is proved. Moreover, as \vn^\^ & V4n^, it follows from (4) that 



lNAV\vNf) 



E 

xeSN„ nr 



V{x)\vNMf 



< 



= m\ 



E 



L 



\vn. 



Hence (6). The estimate (7) is proved in [6]. To prove (8), we notice that 



iiNAv)niNAv)) 



= E 
= E 



{iNAv)){xni^(v)){x) 



\V{x)\' 
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The bound (9) is a straightforward consequence of (7): 

||n2jv.(/iv,(V))||Hi < \\InAV)\\hs < \\V\\hs + 11(1 - In,){V)\\h^ < (1 + Cs,s)\\V\\hs. 



Now, we notice that 



ken* I \k\<^Nc 



^ ( Yl Vk+N.K] efe. (13) 



Prom (13), we obtain 



\U2nA^-In,){V)\\^. 



< 



ken' I |fe| 



E (i + i^i'r 

ken* I |fe|<^-'Vc 

E (i+ifci'r 



E '^k+NgK 

Ken'Xio} 



max 



ifeK*\{o} 



On the one hand, 



E (i + i^H^ 

ken* I |/s|<^JVe 



327r /47r 



N^^oo 2s + 2, \ L 



2s 



and on the other hand, we have for each k e TZ* such that |fc| < ^N^, 



E ^k+N,K 

Ken*\{o} 



< 



c V 

Ken''\{'S) I ^ s I 



< CCol—] N- 



where 



Co = max 



2w J 

y ■ 

yeK3 1 |j/|<i/2 \y - K\"^ 



The estimate (10) then easily follows. Let us finally prove (11). Using (3) and 
(4), we have 



= Ci 



E \ jv 



(l + |^;iVe(x)|^) 

^(i + Kj4) = c^(|r| + ||^^,j|i.J. 



6 



This completes the proof of Lemma 2.1. □ 



3 Planewave approximation of the TFW model 

In the TFW model, as well as in any orbital-free model, the ground state elec- 
tronic density of the system is obtained by minimizing an explicit functional of 
the density. Denoting by J\f the number of electrons in the simulation cell and 

by 



9^AA = |p>0|ypei?^(r), 



the set of admissible densities, the TFW problem reads 
where 

C^fii is a positive real number (Cw = 1, 1/5 or 1/9 depending on the context [8]), 
and Ctf is the Thomas-Fermi constant: Ctf = -^(Stt^)^/''. The last term of 
the TFW energy models the periodic Coulomb energy; for p and p' in H^^{T), 

fee-R.*\{0} 

We finally make the assumption that is a 7?.-pcriodic potential such that 

3m > 3, C > s.t. V/c G U*, < C\k\-"' . (15) 

Note that this implies that 1/'°" is in H"'-^/^^^{T) for all e > 0, hence in C^(r) 
since m — 3/2 — e > 3/2 for e small enough. It is convenient to reformulate the 
TFW model in terms of u = y^. It can be easily seen that 



tTFW 



where 



^TFW(^) = ^ I Vt;|2 + Ctf ^ \v\'^/^ + ^ + \dt{\v\^ \v\''). 

Let F{t) = CTFi^'^ and f{t) = F'{t) = ICtf*^/^- The function F is in 
Ci([0, -^oo)) n C°°((0, 4-oo)), is strictly convex on [0, -hoc), and for all (^1,^2) € 



\f{tl)t2-f{tl)t2-2r{t\)ti{t2-tr)\ < gcTFmax(^l/^^^/^)|^2-^l^ (17) 
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The first and second derivatives of E"^^^ are respectively given by 



where we have denoted by Hj^^ the TFW Hamiltonian associated with the 
density p 

^TFW ^ ^ j^^^ ^ yioa ^ yCoulomb^ 

where 

fceTC*\{0} 

is the 7?.-pcriodic Coulomb potential generated by the 7?.-periodic charge distri- 
bution p. Recall that v^"^^""^^ can also be defined as the unique solution in 
H^^iT) to 



r 



.^y-Coulomb = ( p - |r|-l / p 
I^Coulomb Q 



/ 

Let us recall (see [12] and the proof of Lemma 2 in [3]) that 

• (14) has a unique minimizer p^, and that the minimizers of (16) are u and 
—u, where u = 

• w is in iJ^+^/^"'(r) for each e > (hence in C|(r) since m+l/2-e > 7/2 
for e small enough); 

• u > on M^. 

• u satisfies the Euler equation 

H^^iu) = -^Au + Qctfw'/' + V'°^ + u = Xu 

for some A € K, (the ground state eigenvalue of Hj^^ , that is non- 
degenerate). 

The planewave discretization of the TFW model is obtained by choosing 



1. an energy cut-off > or, equivalently, a finite dimensional Fourier 
space Vn^, the integer being related to E^. through the relation N^. := 
[x/2^L/27r]; 

2. a cartesian grid with step size L/Ng where Ng e N* is such that 
Ng > AN^ + 1, 
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and by considering the finite dimensional minimization problem 

^ (18) 



where 



1 



+ -Dr{\vNA',\vNj% 

denoting the interpolation operator introduced in the previous section. The 
Euler equation associated with (18) can be written as a nonlinear eigenvalue 
problem 

yvN, e Vat,, (CHj^™;^|2 - >^n,,n,)un,,n,,vnJh-\h^^ = 0' 
where we have denoted by 



V 

p 



Coulomb 



the pseudospectral TFW Hamiltonian associated with the density p, and by 
^Nc,Ng the Lagrange multiplier of the constraint /p [uat^P = Af. We therefore 
have 



- — AUN^^Ng+ilN. 



Under the condition that Ng > ANc + 1, we have for all (j) e C^(r), 

v(fc,0 e7^* x7^* s.t. |fc|,|;| < ^iVe, J^iNM4^i = ^k-J, 

so that, n^^^^ is defined on V^^ by the Fourier matrix 

^r (7ijf „ 5 FFT,Ar„ — - — FFT,JV„ 

FFT.ATg 

\k-i? <'-^">- 

where, by convention, the last term of the right hand side is equal to zero for 
k = l. 



We also introduce the variational approximation of (16) 



(19) 
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Any minimizer un^ to (19) satisfies the elliptic equation 



for some Aw e 



(20) 



The main result of this section is an extension of results previously obtained by 
A. Zhou [161. 



Theorem 3.1 For each S N, we denote by um^ a minimizer to (19) such 
that (■UAr<,'^)i| — /'^'^ A^c S N and Ng > AN^ + 1, we denote by 
UNc,Ng a minimizer to (18) such that {uN^^_Ng,u)]^2^ > 0. Then for Nc large 
enough, un^ o,nd UNc,Ng are unique, and the following estimates hold true 





< 




(21) 


lAiVe - A| 


< 


^^^_(2m-l-e). 


(22) 


^ T-TFW T-TFW 


< 




(23) 




< 


(0^jY3/2+(s-l)+jY-m. 


(24) 


1 AjVe, iVg — AjVel 


< 


CN^/'^Ng"'; 


(25) 


1 T-TFW tTFWi 


< 


CN^/'^Ng"', 


(26) 



for all —m + 3/2 < s < m + 1/2 and e > 0, and for some constants 7 > 0, 
Cs,e > 0, Ce > 0, C > and Cg > independent of Nc and Ng. 



Remark 1 More complex orbital-free m,odels have been proposed in the recent 
years [15], which are used to perform multimillion atom DFT calculations. Some 
of these models however are not well posed (the energy functional is not bounded 
from below [2]), and the others are not well understood from a mathematical 
point of view. For these reasons, we will not deal with those models in this 
article. 



3.1 A priori estimates for the variational approximation. 

In this section, we prove the first part of Theorem 3.1, related to the variational 
approximation (19). The estimates (21), (22) and (23) originate from arguments 
already introduced in [3] . For brevity, we only recall the main steps of the proof 
and leave the details to the reader. 
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The difference between (16) and the problem dealt with in [3] is the presence of 
the Coulomb term Dr(|t^p, l^'P), for which the following estimates are available: 

(27) 
(28) 



Q<Dr{p,p) 


< 


CllPlIk, 


for all p e i#(r), 


\D-c{uv,u'w)\ 


< 


C\\v\Ui\ 


|u;||i|, for all («,«;) e (i|(^))^ 


\Dt{p,vw)\ 


< 


C\\p\\li\ 


l^llL|lklli,|, for all {p,v,w) e (i#(r))^ 


II V II i~ 


< 


c||p|U|, 


for all p e i#(r), 


||T/Coulomb|| 

II ^> llj?;+= 


< 


CllPlkj: 


, for all p e H'^{T). 



(30) 
(31) 

Here and in the sequel, C denotes a non-negative constant which may depend 
on r, and M, but not on the discretization parameters. 

Using (27), (28) and the fact that /' > on (0, +oo), we can then show (see the 
proof of Lemma 1 in [3]) that there exist /? > 0, 7 > and M > such that for 

all V G if^(r), 



< ((Wjo™ - \)v,v) < M\\v\\l, (32) 

^||^;||k < ((^^^^"(«) - 2X)v,v) < M\\v\\l. (33) 

and for all v G -ff#(r) such that ||i'||l2^ = AA^/^ and {v,u)i,2^ > 0, 

7lk - ^llk < {('^7^ - ^)(« - - (34) 

Remarking that 

E^™(n^J-i;^™(n) = ((7^J<fw_^^(„^^_^)^(„^^_^)^^_^^^^ 

+ ^Dr{\uNf - \u\\\uNf - \u\'') 
+ [ Fi\u^f ) - Fi\u\') - fiH'KM' - (35) 



and using (34), the positivity of the bilinear form Dr, and the convexity of the 
function F, we obtain that 

jTFW _ jTFW ^ £;TFW(^^j _ JjTFW(^) > ^||^^^ _ ^||2^^_ 

For each A^c G N, mat = A/'^'''^njv„M/||nAr^u||^2 satisfies {un^tU)]^2 > and 

# # 

II^WcIIl^ = A/"^/^, and the sequence (MAr^)^^^^^ converges to u in iJT'^^^^~'^(r) 
for each e > 0. As the functional E'^^^ is continuous on -ff^(r), we have 
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Hence, (wjVe)jVceN converges to u in H^{r), and we also have 



Cw 



, , \Vu^f+ f{\uMf)\ur,f+ VnuNf + Dr{\uNf,\uMf) 



Nc->-oo 



w 



A. 



As f{\uN^'')uN^ + V^'°"ujv. + ^S«j§"''wAfc is bounded in L|(r), uniformly in 
iVc, we deduce from (20) that the sequence (wjVc)jVeeN is bounded in if^(r), 
hence in L°°(r). Now 



^-Coulomb^ 



T/Coulomb^, T/Cc 
+ (1 - n;vj {f{\u\'')u + + 

-Ajv„(MArc - u) - (Ajv„ - A)m . 

Observing that the right-hand side goes to zero in L'^iV) when Nc goes to 
infinity, we obtain that ((iiv,,)A',,eN converges to u in H'^iY), and therefore in 

C°'^/^(r). In addit ion, we know from Harnack inequality [10] that u > in M^. 
Consequently, for Nc large enough, the function un^ (which is continuous and 

7?.-pcriodic) is bounded away from 0, uniformly in N,,- As / G C°°(0, +oo), one 
can see by a simple bootstrap argument that the convergence of {uNS)N^&i to u 



also holds in H 
from (35), remarking that 



m+l/2-e 
# 



(r) for each e > 0. The upper bound in (23) is obtained 



< 



F{\ui,f)-F{\u\^)-f{\u\^){\ui,f-\u\') 



< yCTF^ma^(|«wj4/M«|4/3)|ti;Ve - 



35. 



< ^Ctf (max ||ujvJ|l~ 



4/3 



and that 



< Dt{\unX -H^AunX -\U\')<C\\\UN, 



The uniqueness of ujVe for large enough can then be checked as follows. First, 
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(ujVe,AjVe) satisfies the variational equation 



VwiVe e Vat,, - AjvJujVe,WivJH-\j?^ 



0. 



cs 



Therefore Xn^ is the variational approximation in Vn^ of some eigenvalue of 
?^TFW _ As (TiwJw.GN converges to u in L°°(r), T^^fw^ _ ^tfw ^onverg 

to in operator norm. Consequently, the n^'' eigenvalue of converges 

to the n^^ eigenvalue of Hj^^ when Nc goes to infinity, the convergence being 

uniform in n. Together with the fact that the sequence {XN^N^eJi converges to 

A, the non-degenerate ground state eigenvalue of T-C^J^ , this implies that for 
Nc large enough, A^v^ is the ground state eigenvalue of in Vm^ and for 

all vn^_ e Vn, such that ||wivJ|L| = A/"^/^ and (wat.jWatJl^ > 0, 

+ lF{\vMf) - Fi\uMf) - fi\uMf){\vNf - \uNf) 

> iC^flJ^l^ - AArJ(^^JV, - unJ, {vn, - un,))h-\hI 

> 2\\'>^N,-UN^\Hl■ (36) 

It easily follows that for Nc large enough, (19) has a unique minimizer un^ such 
that («jVe,w)i,2 > 0. 

Let us now establish the rates of convergence of |Ajv<, — A| and ||ujVc — '^^lli?^- 



First, 



(37) 



with 

,„ _ /(["jV^n-Zd^n , |2 , ^Coulomb/ 

As Un^ is bounded away from and / G C°°((0, +oo)), the function wn^ is 
uniformly bounded in least for Nc large enough). We therefore 

obtain that for all < r < m — 3/2, there exists a constant Cr & K+ such that 
for all Nc large enough, 

|Ajv. - A| < a (\\uN^_ - u\\li + \\uN^ - • (38) 

In order to evaluate the H^-norm of the error {u^^ ^ ^*) , we first notice that 

yvN, e Vat,, ||wjv, - < ||wjv, - vnJ\h^ + \\vn, - u\\h^, (39) 
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and that 

# # ' # 

= - 2X)iuN^ - u), {UN^ - vnJ)h-\h^^ 

- 2X){u - v^J, {ur,, - t^ArJ)^-!,^^) • (40) 

For all zpf, e Vn^, 

((i?T™"(u) - 2X){u^^ - «),^ivJH-,Hi 

= -2j^[f{\uNA^)uN^-fi\u\^)u-2f'i\u\^)\uf{UN^-u)]zN^ 

-2Dr{{uN, - u){uN, + u), (ujv, - u)znJ - 2Dr{{uN, - uf,uzN,) 
+2{Xn,-X) I UN^ZN,. (41) 



On the other hand, we have for all vn^ S Vn^ such that Hwat^Hl^ = M^I'^, 



# 



Using (17), (29), (38) with r = and the above equality, we therefore obtain 
for all fjVe S Vn^ such that ||fjvJ|L| = -A/"^^^, 

((^^^^"(m) - 2X){ur,^ - u), (UN^ - vnJ)h-\h^^ 



+ (ll^iVe - wllff^ + Wun^-uWli) \\un, - VNj\l2^y (42) 
Therefore, for Nc large enough, we have for all vn^ £ Vn^ such that ||t^Ar<.||z,2 = 

Together with (39), this shows that there exists iV G N and C G M+ such that 
for all iVc > N, 

yvN^ G Vn^ s.t. \\vnJ\lI II^Wc - wllff^ < C\\vN^-u\\Hy 

By a classical argument (see e.g. the proof of Theorem 1 in [3]), we deduce 
from (2) and the above inequality that 

\\uN^-u\U<C min \\vn^-u\\hi <C^,,N-^^-'/^-'\ (43) 
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for some constant Ci,r independent of Nc- This completes the proof of the 
estimate in the if^-norm. We proceed with the analysis of the L^-norm. 



For w e i^(r), we denote by ipw the unique solution to the adjoint problem 
f find ipyj e U"*" such that 



(44) 



where 



The function ijjyj is solution to the elliptic equation 



from which we deduce that if w € H^{r) for some < r < m — 3/2, then 
e H^/^T) and 



< C'r||w||ijr , (45) 

# ^ 



for some constant Cr independent of w. Let u*j^ be the orthogonal projection, 
for the inner product, of un^ on the affine space G -^#(r) | J^uv = J\f^. 
One has 
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from which we infer that 



= ^("Afc - -U) + ^\\un, - All (^^- J^'^N,^ 

= {UN, - u,u%^ - u)h-\hI + 4^ll^^c - All 



= {{E {u) -2\)ijuM,-u,UN^-u)jj-i^^^^ + ^\\UN,-A 



= {{E {u) -2X){un,-u),iPum,-u}h-\hI + ^^No-AlI 



1 

. . rW'^Nr — Ar^ 



= {u) - 2X){UN, - «), V^„„,-«)^^^i,^^ + ^ll«A^. - Ai^ 

+jf\\uN, - All j^f i^'^)^^^»N,-u + Dr{u^ ,u^uMc-u) 
For all tpNc ^ VaTc) it therefore holds 
ll"iVe-w|||2 = {{E'^^^"{u) -2X){un^-u),iPun,-u-->PnJh-^ 



# ' # 

1 



(u) - 2X){uN^ - u), V^J^-i „i + -r.lli.iv. - All 



(46) 



Using (17), (29), (38) with r = and (41), we obtain that for all ipN, e Vn, Hu-^, 



{{E^'^^ {u) - 2X)iuN^ - u), ^nJh-\hI 



< C{\\UN^-All 



+ \\un, -AlI [\\un, -AhI + W^N, -Ali) llV'iVcllff 



(47) 



Let us denote by 11^^ ^^^j^ the orthogonal projector on Vat. D u-^ for the 
inner product and by = 11^^ j^^^-0«iVe-M- Noticing that 



U%\\Hi < Uu^^^uWhI < r'M\\uN^^ - U\\ 



we obtain from (33), (46) and (47) that there exists A*' e N and C e M+ such 
that for all Nr > N, 



hN,-Ali < c l\\uN,-ALi\\uN,-AHi+hN,-AHi\\i^u^,-u-i^°Nj\Hi 
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Lastly, for all v € u-^ and all A^c e N* 

ll--n^„^n„..|Ui < (^ + 2^LWJ,» ) ll^-n-c-ll^i' (48) 

so that, in view of (2) and (45) 



< CN-^UN^-u 
Therefore, 

|hAr„-u||L| < C(^\\UN^-U\\l^^+N-^\\UN^-U\\H1^ 

By means of the inverse inequality 



/27r 



(r-s) 



yvN^eVr,^, |KJk;<(^-^J iV/-^IK.Ikj, (49) 

which holds true for all s < r and all A^c > 1, we obtain that 

\\UN^ - u\\h^ < Cs,,N-^"'-'+^/^-'^ for all < s < m + 1/2. (50) 

To complete the first part of the proof of Theorem 3.1, we still have to compute 
the H^'^-novm of the error (ujVe — u) for < r < m — 3/2. Let w G H^{T). 
Proceeding as above we obtain 

w{uN,-u) = ((i;'^™"(u)-2A)(uAr^ -u),n^^ 
Jr = # ' # 

+{{E'^^'^"{u) - 2X){UN^ - u), V. - n^v^^nu^i'^v)H-\Hl 



l^uw. (51) 

Combining (33), (45), (47), (48), (50) and (51), we obtain that there exists a 
constant C G 1R+ such that for all A^c large enough and all w G H^{r), 

J wiUN^ -U) < C (ll^AT, - U\\l,^ + iV,-('-+l)| 
< C_.,,7V-(™+'^+V2-^)||^||^3.. 



UN, - \\W\\H, 



# 
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Therefore 

\WN^-u\\H-r= sup -^^^ <C_,,,7V-('"+^+V2-), (52) 

* »eHj(r)\{o} \\w\\h^ 

for some constant C_r,e € IR+ independent of Nc- Using (38), (43) and (52), we 
end up with 

lAjv, - A| < C,iV-(2™-i-^). 



3.2 A priori estimates for the full discretization. 

Let us now turn to the pseudospectral approximation (18) of (16). First, we 
notice that 

< £;TFW(_^l/2|r|-l/2) 

< CTFAA-'^/^|rr2/3 + ||r°"||L-AA, 

from which we infer that ujv^jVg is uniformly bounded in H^(T). We then see 
that 



2 



+-Dr(|wjVc,jv„P, l^jVe.JvJ^) • 



Using (6), (11) and (27), wc obtain that Xn^.n^ ^^l^o is uniformly boimdcd. Now, 

+2C^^njv„ (V|^°"|°™pMjVe,JVg) - 2C^^XN„NgUN„Ng, (53) 

and we deduce from (4), (6) and (8) that 

\\^N,i^N,{f{\uN,,Nj^)UN,,N,))\\^2^ < U (^^iV, (/(|wjV„,JV, T))) \un,,N, 



-CTF||wjVc,iVg llii ||wjVe,JVg 11^4^, 
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and that 



lin^v. {In, {V'"^ur,^,N,)) hi < \\^2n^ {In, {V'°^un^,n,)) 



# 

1/2 

ion|2i . |2\ 



Besides, using (30), 

As UNc,Ng is uniformly bounded in H^{r), and therefore in -^^(r), we get 

\Wn,,nJhI = (ll"Wo,wJli^ + II^WWc,AfJli|j 

< c(i + iK„ivjii/i) 

< c(l + |K„;vJ|^|). 

Therefore un^,n^ is uniformly bounded in H^{r), hence in _L°°(M'^). 

Returning to (53) and using (9), (15), and a bootstrap argument, we conclude 
that UAf^^jVg is in fact uniformly bounded in H'^'^~^'^{T). 

Next, using (36), 

= i;^™(«iv.,^j-i?^r(«ivj 

+ l^{l-I^^)iF{\uM^,^f)-F{\u^f)) 

< l^{{i-ir,j{vm\u^^,r,f-\uNf) 

+ 1 (i-i„)(F(K^„ n-i^(Kj2)). 



Let g{t, t') = t'-t^* ^ ■ large enough, un^ is uniformly bounded away 

from zero; besides, both un^ and UN^,Ng are uniformly bomidcd in h'^'^'^'^ (T) . 
Therefore, ff(wjVe, wjVe,JVg) is uniformly bounded in H^^~^'^{T). This implies that 
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the Fourier cocfRcicnts of g{uN^,UN^,N ) go to zero faster that \k\ which in 
turn imphes, using (5) and (10), that 

^(1 - X^J(F(|«;V.,ivj') - F(|«ivj')) 

= "^^s) (5("JV,,Miv,,JvJ) {un^^n, - UnJ 

< CN^J^N-y^ \\uN^,r,, - . (54) 

On the other hand, 

{{l-I^^){V)){\u^^^,r,f - \uNf) 
< \\^2nA{^ -1-n,){V))\\l2J\un,,n, +un,\\loo\\un,,n, - unJlI 



Therefore, 



(55) 



We then deduce from (55) and the inverse inequahty (49) that {uN,^,Ng)Na,Ng>iN^+i 
converges to u in -ff|t(r), and therefore in L°°(]R^). It follows that for Nc large 
enough, ujVo.Wg is bounded away from zero, which, together with (53), implies 

that {uM^,Ng)Nc,Ng>ANc+i is bounded in H^'^^^'^~'^{V). The estimates (54) and 
(55) can therefore be improved, yielding 



L 



T 

and 



{l-lNMF{\uM^,j,f)-F{\uNf)) 



<CNl''N-^"^^y'-nuN.,^g-u^^ 



We deduce (24) from the inverse inequality (49). For Nc large enough, UN^^Ng is 

bounded away from zero, so that fi\uN^,Ng P) is uniformly bounded in H^~^^^'^~'^{T). 

Therefore, the k*"^ Fourier coefficient of (y°" + .f{\uN^,Ng\'^)) is bounded by 
C|fc|~"* where the constant C does not depend on Nc and Ng. Using the equal- 
ity 



A 



-A 



= J\f- 



iC^IuJ^l^ - AjvJ(wjVe,JVg - UN J, iUNc,Ng - UNc))h-\h\ 



- y^(l-I;vJ(V'- + /(|«;Ve,ivj'))K. 
U{\uNc,Ng\')-f{\uNmUNc,Ng? 



Na 
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(24) and (29), we obtain (25). A similar calculation leads to (26). 
Lastly, we have for all vn^ G Vnc^ 



£^^r(^;vj - i?^7(^.^.,^.J (56) 

= ( C^ul^N, - Wc, ) (^'A'c - WAT,, AT J , (VN, -UN,,N,))H-\Hi 

+ ^-Dr(|vjvJ^ - |MJV„,JvJ^ |vjvj^ - |mjVc,JvJ^) 

3 



+ E iw) (ni^iv.(^)|')-i^(|nA,,(a;)n-/(|n;v.(a;)|^)(|^;A,,(a;)|2-|«;v.(a;)|2)) 



As UN^,Ng converges to u in i?|t(r), the operator "^j^^^p — 'M.J^^ converges 

to zero in operator norm. Reasoning as in the proof of the uniqueness of Un^, 

we obtain that for A^c large enough and Ng > iNc + 1, we have for all Vn^ G Vn^ 

such that ||fjv<.||L2 = A/"^^^ and (vnciUnJl^ > 0, 
# # 

Thus the uniqueness of UNc,Ng for Nc large enough. 

4 Planewave approximation of the Kohn-Sham 
LDA model 

The periodic Kohn-Sham LDA model with norm- conserving pseudopotentials [14] 
leads to the constrained optimization problem 

7^^ = inf {E^^i'i?), $ e A1} (58) 

where 

= 1$ = (01, • • • , M'^ e (77^(r))^ I ^ <f>,<t>j = 5,j^ , 

Af being the number of valence electron pairs in the simulation cell, and where 

£;^S($) = V / |v</.i|'+ / p$viocai + 2V(0i|Ki|0i) + J(p*) + i?^°^(p$). 

(59) 
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The density /?$ associated with <I>, the Coulomb energy J{p<i,) and the LDA 
exchange-correlation energy E^^^{p,j,) are respectively defined as 

AT 

J{p^) = ^£»r(p$,P*) = 27r ^ \k\-^\(f^)k\^ 

feeK*\{o} 

pLDA/ ^ _ / „LDA 



i^xc (P*) = l^e^!f''{Pc{x) + p^{x))dx, 



where pc > is the nonlinear core correction and where e^^^{p) is approxi- 
mation of the exchange-correlation energy per unit volume in a uniform electron 
gas with charge density p [8]. 

The local and nonlocal contributions to the pseudopotential model the interac- 
tions between valence electrons on the one hand, and nuclei and core electrons 
on the other hand. Troullier-Martins pseudopotentials [14] constitute a popu- 
lar class of pseudopotentials for which the Fourier coefficients (Viocai)fe decay as 
with m = 5. The nonlocal contribution is defined by 



M 



where the functions Xj regular enough functions of L|^(r). In all what 
follows, we will assume that 

3m > 3, C > s.t. Vfc e TZ*, |(lWi)fe| < C|fcr™ (60) 

and that 

vi<i<M, Ve>o, xi e F^-'/'-'(r). (61) 

The function p n- e^^^(p) does not have a simple analytical expression. Al- 
though this function is of class C°° on the open set (0, -)-oo), DFT simulation 
softwares make use of approximate functions which are C°° on (0, p*)U(/9*, -|-oo) 
but only C^-^ in the neighborhood of the density p^ := 3/(47r) (atomic units) [8]. 
In order not to deteriorate the convergence rate of the pseudospectral approx- 
imation, it is better to ressort to more regular approximations of the function 
gLDA ^ggg jgj^ .^jjj assume here that 

the function p ^ e^f ^(p) is in Ci([0, -|-oo)) n C['"l((o, -^oo)), (62) 

J LDA 

e^?» = 0, ^(0) = 0, (63) 

(where [m] denotes the integer part of [to]) and that there exists < a < 1 and 
C e K+ such that 



Vp G M+ \ {0} , 



2„LDA 



XC 



-ip) 



LDA 



rfp3 



-ip) 



<C(l + p"-^). (64) 
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Note that the Xa exchange-correlation functional {e^"{p) = — Cxp'*^'', where 
Cx > is a given constant) satisfies the assumptions (62)-(64) with a = 1/3. 
Let us also remark that (62) and (64) imply that 

e^°^ e C^'" ( [0, i] ) for each L > 0, (65) 

a property we will make use of below. Lastly, we assume for simplicity that 

Pc€C^(r). (66) 

It is easy to prove that under assumptions (60)-(66), (58) has a minimizer $° = 
{<j)i,--- ,<j>%-y with density = p^o. The regularity assumptions on Viocah 
on e^^^ and on the functions Xj allow to state that the minimizer is in 
[H%{T)]^, and even in [if™+^/^"'(r)]^ for any e > 0, if at least one of the 
following conditions is satisfied: e^^^^^ G C[™1([0, +oo)) or + > in M? . The 
former condition is not satisfied for usual LDA exchange-correlation functional. 
On the other hand, it is satisfied for the Hartree (also called reduced Hartree- 
Fock) model, for which e}^^ = 0. The latter condition seems to be satisfied in 
practice, but we were not able to establish it rigourously. 

Let us introduce the Kohn-Sham Hamiltonian 



1 / //pLDA 



dp 



h + VpO 



where 



and 



/l = -^A + yioeal + Kl, (67) 
J LDA 

We notice that £;KS'($o^ ^ 4-H^s$o -^^ (_ff-i(r))-^ and thus the Euler equations 
associated with the minimization problem (58) read 

where the Af x J\f matrix = (A^^ ), which is the Lagrange multiplier of the 
matrix constraint Jp (picpj — 6ij, is symmetric. 

In fact, (58) has an infinity of minimizers since any unitary transform of the 
Kohn-Sham orbitals $° is also a minimizer of the Kohn-Sham energy. This is a 
consequence of the following invariance property: 

V$ e A^, V;7 e U{U), U^&M and E^^(C/$) = E^^{^), (69) 
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where U{M) is the group of the real unitary matrices: 

U{M) = {U& M^^^ I U^U = , 

Itv" denoting the identity matrix of rank M ■ This invariance can be exploited 

to diagonalize the matrix of the Lagrange multipliers of tlic^ orthonormality 
constraints (see e.g. [8]), yielding the existence of a minimizer (still denoted 
by with same density p°, such that 

K^^l = 4't>l (70) 

for some e? < < • • • < e^. 



Remark 2 The Kohn-Sham Hamiltonian 'H^o is an unbounded self-adjoint op- 
erator on i^(r), hounded below, with compaet resolvent. Its spectrum therefore 
is purely discrete. More precisely, it is composed of a increasing sequence of 
eigenvalues going to infinity, each of these eigenvalues being of finite multiplic- 
ity. It is not known whether ej, are the lowest eigenvalues (counted with 
their multiplicities) of (Aufbau principle). However, it seems to be most 
often (though not always) the case in practice. On the other hand, the Aufbau 
principle is always satisfied for the extended Kohn-Sham model, for which the 
first order optimality conditions read 



+ 00 

where ep is the Fermi level (see [4] for details). In this article, we focus on the 
standard Kohn-Sham model with integer occupation numbers. We do not need 
to assume that the Aufbau principle is satisfied, hut our analysis requires some 
coercivity assumption on the second order condition at $° (see (73)). 



+ 00 



p°(x)=2^n,|<^°(a;)|2 

= (5,,-, l<i,j<+oo 
ni = l ife°< €F, ni = Oife°> cp, 



For each $ = (^i, • • • , (f>M)'^ € A4, we denote by 

T^M = |(V^i, • • • , iJMf e (i^#(r))-^ I ^ (j^ii^j + ^i<^j = o| 

the tangent space to at and by 

= I* = (v,!,.. . ,i;^f G (iJ^(r))^ I J^Mj = o| . 
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Let us recall (see e.g. Lemma 4 in [13]) that 

where A= {A€ R-^^^ | = -^} is the space of the M xM antisymmetric 
real matrices. 

Since the problem we are considering is a minimization problem, the second 
order condition further states 



'iW e T^oM, 04,0 {W, W) > 0, 



where 



^ i=i -^r 



N 

(71) 



r J2 LDA 

+4E / -^{p. + p'mi<kW (72) 

It follows from the invariance property (69) that 

04,0 (*,*) = for all * e ^i>°. 

This leads us, as in [13], to make the assumption that a^a is positive definite 
on $0'-^ , so that, as in Proposition 1 in [13], a$o is coercive on (for the 

norm). Thus, in all what follows, we assume that there exists a positive 
constant c$o such that 

V^-e^O'-U-, a4,o(*,*) > c$o||*||^i. (73) 

# 

In the linear framework ( J = and E^^^ = in (59)), this condition amounts 
to assuming that there is a gap between the lowest A/""^ and (A/"+l)^* eigenvalues 
of the linear self-adjoint operator h — —^A + Viocai + Vni- 

The planewave approximation of (58) reads 

I^Ij,^ = inf {E^^^i^Nj, $iv. e < n A^} , (74) 

where 



i=l j^;^ 

+J(p$) + ^Xiv,(e^°^(pc + n2iv.P*)). (75) 
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Here Nc is a given positive integer, equal to [^/2Ec L/2tt], denoting the so- 
called cut-off energy, and Ng > ANc + 1 is the number of integration points per 
direction used to evaluate the exchange-correlation contribution. The energy 
E§^{^) is defined for each ^ e M. For <i> gV^ nM, Hsat^ p$ = p$5 so that on 
this set, E^^ differs from E^^ only by the presence of the Fourier interpolation 
operator In^ in the exchange-correlation functional. Let us mention that in 
practice, the terms involving the local and nonlocal components of the pseu- 
dopotential are also computed by some interpolation procedure. However, these 
terms are calculated using spherical harmonics and a very fine one dimensional 
radial grid, so that the resulting integration error is usually much smaller than 
the interpolation error on the exchange-correlation term. Note that, in addi- 
tion, the pscudopotential gives rise to linear contributions that can be computed 
very accurately once and for all (and not at each iteration of the self-consistent 
algorithm). We postpone the analysis of (74) to a forthcoming article [5], and 
focus here on the variational approximation 

= inf {^kS($nJ, $iv. e n M} (76) 

of (58). The unitary invariance of the Kohn-Sham model must be taken into 
account in the derivation of optimal a priori error estimates. One way to take 
this invariance into account is to work with density matrices (see e.g. [4]). An 
alternative is to define for each ^ G M the set 

:= l"^ eM\\\^ -Mr2 = min ||[/* - $11 r2 I , 

and to use the fact that all the local minimizers of (76) are obtained by unitary 
transforms from the local minimizers of 

= inf {E'^^i^Nj, $;v, eV^^nM^"}. (77) 
The main result of this section is the following. 

Theorem 4.1 Assume that (60)-(66) hold. Let he a local minimizer of (58) 
satisfying (73). Then there exists r° > and N° such that for Nc > N°, (77) 
has a unique local minimizer in the set 

{4>^„ e < nAl*" I 11$^^ - $°||^. < r°} . 

// we assume either that e!^^^ € C[™1([0, -|-oo)) or that pc + p° > on T, then 
we have the following estimates: 

\\K^-^°\\h^ < Cs,eN-^'^-'+^'^-'\ (78) 
\4,N.-4\ < C,7V-(2"-i-^), (79) 
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for all —TO + 3/2 < s < to + 1/2 and e > 0, and for some constants 7 > 0, Cg.t > 
0, Ce > and C > 0, where the e°jy^ 's are the eigenvalues of the symmetric 
matrix A%^, the Lagrange multiplier of the matrix constraint /p (pi,N^(pj,Nc = ^ij- 

4.1 Some technical lemmas 

For $ = ,</.^f G (^^^(r))-^ and * = (V'l,--- -V-A^)^ G (^^#(0)-^, we 

denote by M$ $ the J\f x J\f matrix with entries 



-1/2 



r 

The following lemma is useful for the analysis of (77). We recall that if A 
and B are symmetric N x N real matrices, the notation A < B means that 
x'^Ax < x'^Bx for all x e M^. 



Lemma 4.2 

1. Let^ G and* G M. IfM^ii^i^ is invertible, thenUq,^^, = M^^^{M^^^M^^^) 
is the unique minimizer to the problem m.mu^u(^j^'j \\U^ — *IIl|j- • 

2. Let^ eM. Then 

A^* = {(W - Mw,wy^^'^ + W \ W 0< Mw,w < W} 

where l^v denotes the identity matrix of rank J\f. 

3. Let $ = {(f)i, ■ ■ ■ , ^aa)"^ € M. //TVc e N ^5 such that 

dim(span(njv>i, • • • , UN^(f)_sr)) = U, 
then the unique minimizer of the problem min^ &v^nM II^JV, — ^IIl^ is 

7r^^$ = (Mn^^*,n^^*)"'/'n;v,$. (81) 
In addition, ttJ^^ e VjO^ n Al*, 

- < \/2||n;v,$ - (82) 

and for all la,rge enough, 

IK* - Hh^ < ll*IUil|niv.* - *|li^ + ||n^.* - Hny (83) 

4. Let Nc such that dim{VN^) > and $jVe € Vj^ n At. Then 

y^nAl*"^ =[{lN- Mw^^,w^J^'^<^N^ + Wn^ I Wm^ G n (l>]^^, < Mw^^,w^^ < 
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Proof In order to simplify the notation, we set M = For each U e 

WU"^ - ^lll^^ =2J\f-2Tr {MU). 
Any critical point U of the problem 

max Tr (MU) (84) 

satisfies an Eulcr equation of the form At/-^ — M for some symmetric matrix A. 
Besides, Tr (MU) = Tr (A) and A^ = MM'^. Any maximizer U of (84) therefore 
satisfies M = {MM'^Y^'^U^ . Consequently, if M is invertible, the maximizer 
of (84) is unique and reads U^,^^ = M'^ {MM'^)-^/^ . It also follows from the 
definition of the matrix M that * = M$ + 11^ with W e^-^. Thus, 

f7*,4,* = M^(MM^)-i/2^$ + W, 

with W = U^^^W e 

Let us now prove the second statement. Each e {H^{T))-'^ can be written 
as ^1 = + W for some matrix M e R-^^^ and some W e A simple 
calculation leads to 

j^i^iipj = [MM^]ij + [Mw,w]ij- 

Hence ^ = + W e M ii and only if MM^ + Mw,w = W- In addition, 
* G 7W* if and only if * € Al and Uq,,^ = M^{MM^)-^/^ = Ij^, that is to 
say if and only if M is symmetric, < M-^^^ < Ij^ and M = {Ij^ — Mw^wY^^. 

Let (X(u)i</:t<dim(Vjv„) be an orthonormal basis of Vn^ (for the inner product) 
and let C e ]R<i™(^"c)x^ be the matrix with entries 



Note that 



d^m{V^rJ 

For all $jv<: = i4'N^,i, ■ ■ ■ , 4'N^,j\r)'^ G "l^iv^ fl A^, each (/)jVc,i can be expanded as 

dim(VNj,) 

^i-JVe.i = ^ C',*iX;*, (86) 

where the matrix C = [C^,] G x A/" satisfies the constraint C^C = l^A- 

The expansions (85) and (86) can be recast into the more compact forms 

Um^ = C^X and = C'^X, 
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where we have denoted hy X = (xir"" iXdini(V]v ))^' ^ simple calculation 
then leads to 

11$^^ -$|||, =2Ar-2Tr(C^C). (87) 

# 

Reasoning as above, we obtain that the unique solution to the problem 

max Tr (C^C) 

is C = C(C^C)~^/^. Note that the rank of the matrix C is N provided that 
dim(yjv^) is large enough so that the matrix C^C is invertible provided that 
dim(VjVe) is large enough. As a consequence, the unique solution to the problem 
mm^^^lyu^j^ ll*^^ - $11^^ is = (C^C)-^/KFX = {C^C)-^/^Iir,^'^. 

It is then easy to check that C^C = -Mnjv^$,njvc*- Hence (81). Then, for all 
U g M^xAA g^(,ij ^1^^^ ijT^ ^-^j^^ 

||C/<$ - = 2(1 - Tr (t/M^^^^n.,*)). 

and the same argument as above leads to the result that this quantity is min- 
imizcd for U = MA/„%,n.,*(A^n„%,n„,*^n„%,n.,*)-'/' = W- Therefore, 

We also infer from (87) that 

IK^ - ^lli^^ = 27V - 2Tt {iC^Cf/^) = 2Tr - [C'^Cf'^) • 
Besides, an easy calculation leads to 

llHiv,* - $11^2^ = Tr (W - C^C) . 

Using the fact that 

< (w - {c^cf'^) < (w - ic^cy/') (w + (C^C)i/2) = w - c^c, 

we obtain 

\\n^^^-m.^ = 2Tr (l^ - (C^C)^^) < 2Tr (l^ - C^c) = 2||n;Ve$-$|li|. 
Hence (82). We also have 

= ||((Mn„^4.,n„,*)-'/' - 1^)UnM\h^^ + llnw,* - $11^^ 

< !l(Mn,,^*,n„,*)-i/' - wilF||nAr^$||^^ + lin^r^* - $11^^ 
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where || • ||f denotes the Probenius norm. We then notice that 

Consequently, for Nc large enough, 

||(Mn„^$,n„^$)"^/^ - WIIf < ||Mn„ a,-*,n„ *-4.||f < W^n,^ - ^Ili2 . 

# 

Therefore (83) is proved. 

Lastly, the fourth assertion easily follows from the second one. □ 
Lemma 4.3 Let 

K = {W€ {LliT)^ I < Mw,w < W} , 

andS : M^^-^ (the space of the symmetric J\f x J\f real matrices) defined 
by 

S{W) = (Iaa - Mw,wY'^ - W- 

o 

The function S is continuous on K and differentiable on the interior K of K. 
In addition, 

^WeK, ||W||F<||l^||i., (88) 
and for all {Wi,W2,Z) € K x K x (i|(r))^ such that < | and 

\\W2hl ^ ^' 

\\S{W^)-S{W2)\\f < 2(||W^i|U^ + ||W^2||l|)||W^i-W^2|U|, (89) 
\\{S'{W,)-S'{W2))-Zy < 3\\Wi-W2hi\\Zhi, (90) 
\\{S"{W^){Z,Z)\W < S\\Z\\l. (91) 

Proof Diagonalizing Miv,w and using the properties of the function i i-)- (1 — 
^^1/2 _ I t^q ggg ^iia^ ^ ig continuous on K and differentiable on K, and that 

IIWIIf < \\MwMf < \\w\\li- 

Hence (88). As 

S{W) + ^S{Wf = -^Mw,w, 

o 

we have for allW G K, 

S'{W) ■Z+^ [S{W){S'{W) ■ Z)) + {S'{W) ■ Z))S{W)] 
= -^[Mw,z + Mz,w]- 
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Denoting hy A = S'{W) ■ Z, we deduce from the above equality that 
\\A\\l+Tr{A^S{W))<\\AMMwMW<\\A\\^\m^^^^^ 

As \Tr{A^SiW))\ < \\A\\US{W)h < \\A\\l\\SiW)y < \\A\\UW\\l, we fi- 

# 

nally obtain the inequality 

\\AMl-\\W\\l.J<\\Wh.jZh.^, (92) 

which straightforwardly leads to (89) under the conditions ||V[^i||z,2 < ^ and 
II 1^211^1 < I Lastly, 

{S'{W2) - S'iWi)) ■Z + ^ [S{W2){{S'{W:,) - S'{Wi)) ■ Z) + {{S'iW^) - S'{Wi)) ■ Z)S{W2)] 

+ 1 [{S'iWi) ■ Z){S{W2) - S{Wi)) + {S{W2) - S{Wi)){S'{Wi) -2)] = -^ [Mw,-w„z + Mz,w,-wA , 
so that still under the conditions < 5 and ||W2||i,^ < 5, 

^||(5'(w^2) - s'{w,)) ■ z||f < ^||w^2 - i^ilU^ll^iu^. 

Hence (90). Lastly, taking W2 = Wi + tZ in (90) and letting t go to zero, we 
obtain (91). □ 



Lemma 4.4 Let be a local minimizer of (58) satisfying (73). Then a^a 
defines a continuous bilinear form on (iJ^(r))-^ x {H^(r))^ , and there exists 
N* such that for all > N*, 

lK$o_$0||^, <1, (93) 
a<,o«$o - - $0) > ^|K$° - ^Yh^ , (94) 

VWe[<$T, a^o{W,W)>^-f\\W\\l.^. (95) 

In the sequel, we denote by C$o the continuity constant of a$o , i.e. 

V(*,*')e((^^^(^))^)^ \a^o{^,^')\<c^o\mHiW\\Hy (96) 

Proof Estimate (93) immediately results from the closeness of tt^ $0 to 
Using the fact that n-^^^^ S At*" (see Lemma 4.2, point 3), we get 

7r^$o_$o =5(W^)$o^p^ (97) 
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with W e [$'']"^, from which we derive, using (88), that 
04,0 (77^^^$° - $°,7r^^$° - $°) = a^o{W,W) + 2a^o{W,S{W)^^) + a^o{S{W)^°,S{W)^°) 



> c^o||W||^. -2C^o||T^||^. -C4.o||W^||td|$' 

# 9P # # 

> (c^o-2C7^o||H^|U^||$"||^. -C^o||W^||i^||$0||^,)||W^||^, 



As by (82), ||7r;)^$° — '^'^IIl^ goes to zero when Nc goes to infinity, so does 
\\W\\l2 . Using again (88), we deduce from (97) that 11 11^1 ~ Ikw ^° - 
Hence (94). 

Finally, for each W G [tt^^^O]-^, W^* = - Mvi.,4,o$o belongs to [$°]-^. Re- 
marking that Mvy,4>o = My^i^o_T^Mi^o, we derive 

||Miy,<i.o||F < \\Mw,^o_^M^4y < e(iVc)||VF||L| 

where £{Nc) = ||$° - 7r5^$"||^2^ when iVc goes to infinity. Therefore, 
As 

0*0 {W,W) = a*o ( VF* , VF* ) + 2a$o {W,W -W*) + a^o{W -W* ,W -W*), 

we obtain 

a*o(W^,lV) > c$o||M/*||^^-2C$o||VF||^^||W^- 1^*11^^ -C$o||W^-W^*||^^ 

> (c$o - 2C*0£(Are)||$°||ff^ - C^oe{N,f\\^^\\H^^ WWfjj,^. 

Hence (95) for large enough. □ 



Lemma 4.5 There exists C > such that 



1. for all (Ti, T2, T3) G ((^^(r))^' ' 



2. for all Ti e (/fit(r) n L^iT))^ and (T2, T3) G ((i?^(r))^)', 

(i?KS"(<i," + Ti) - i?KS"($o^^ (T^^^ T3)| < C (1 + IITilli"") ||Ti||2^ IIT^Il^. ||T3||^i 
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Proof Let us denote by 

Splitting r$o (Ti, T2, T3) in its Coulomb and exchange-correlation contribu- 
tions, we obtain 

r$o(Ti, T2, Ts) = rg?"(Ti, T2, T3) + r-o(Ti, T2, T3), 



with 

„Coulomb 



rgo°"(Ti,T2,T3) = 16^ (i?r((/'S,i,t;2jt^3j) + £>r(0S,i.^iJ^^3,,) + i5r(</'^3,i,t^i,jt^2j 

1.3 = 1 

+16 ^ -Dr(wi,iU2,i,wi,jU3,j) + 8 ^ DT{vl^i,V2,jV3^i), 



and 



r|'6(Ti, T2, T3) = rl'o'ir,, T2, T3) + rJ^o''(Ti, T2, T3), 



where 



' $0 



^(Ti,T2,T3) = 4^(^ 
rJo''(Ti,T2) = 16^ 



dp 



i2„LDA 
xc 



)A / \ / \ -, 



Using (29), wc obtain that there exists a constant C > 0, such that for all 

(Ti,T2,T3)e ((iJ4(r))^)', 

|^Couiomb(Y^,T2,T3)| < C (lITilU^ + IITilli.J IIT2IIH1 IIT3IU. . (98) 



Using (64), we get 

(Pc + P*o+Ti ) 



dp 



from which we infer 



(Pc +P$o) 



< C 



a/2 , 

Pr\ +Pri 



|r^^o'^(Ti,T2,T3)| < (pxf P^'IV^/'. 



(99) 
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Introducing the function 



we can rewrite r^d [i-i, 1 2, is) as 



r^o'(Ti,T2,T3) 



16 



= 16 



Ir Jo 



j2 LDA 
"* '^xc 

1 rd2eLDA 



dp 
dp- 



(Pc+P$(l)) {^J)i{^)v2^ 

\i=l / \i=l 



^2gLDA 



^3gLDA 



+2^p-(Pc + P$(t)) \^^<i)i{t)Vi,ij \^^<i)i{t)V2,ij \^^<j)i{t)V:i, 

Thus, using again (64), we obtain 



r^f (Ti,T2,T3)| < Cj j {l + {Po+Pmr'')PHt)Pr:pr:P'4. 



a-lx 1/2 1/2 1/2 1/2 



dt 



- ^Ll!^^^ Pm^pl%py>\^>r^s dt- 



Now, 



\i=i i=i i=i 



N \"-i/2 

■2 1 



•^0 \i=l Z^i=i^l,i 



a-1/2 



dt. 



< 2«-i/2 
Therefore, 



,^E£i0K 



2a-l 



-1/2 



^i=l 



dt < 



a-1/2 



|r-'^(Ti,T2,T3)| < C|(p:?:"("'^/^)+pt,Jpx'>t' 
Gathering (98), (99) and (100), we obtain the desired estimates. 



1/2 1/2 



(100) 
□ 
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Lemma 4.6 Let be a local minimizer of (58) satisfying (73). Then there 
exists C > such that for all E Ai , 

E^^{^) = E^^{^°) + 2a4,o(* - * - $°) + - $°), (101) 

with 



Proof Using the fact that the first order optimality condition (70) also reads 
[^KS'(^O)]. ^ 47^KS^o ^ 4g0^o in if-i(r), we have for all * G 7W, 

£;^^(*) = i;^^($°) + (i;^^'($°),*-$°)^-i^i +;^£;^^"($°)(*-$°,*-$°) 

+ / (£;KS"(^0^^(^_^0))_^KS"(^0))(^_^0^^_^0)(^_g-)^g 
AT 

= E''^{<^')+4y^e^ / 0O(v>i-0O) + ii?KS"(cj>O)(^_$o^^_^O) 
i=i -^r ^ 

+ / (£;KS"($o + s(* - $0)) - E^^"{^°)){^ - * - $°) (1 - s) ds 
Jo 

= i;KS($0) _ 2 V e° / (Vi - .^°)' + ^£;^s"($°)(* - $0, * - $0) 

+ / (£;KS"((^O^g(^_^0))_^KS"(^0))(^_^0^^_^0)(^_^-)^g 

= E^^{^'^) + 2a$o(* - * - $°) + i?(^' - $°), 

where 

i?(T)= / (£;KS"($o + sT)-£;^s"(^0)^(y^^^^^_^^^^_ 
Ja 

The estimate (102) then straightforwardly follows from Lemma 4.5. □ 



4.2 Existence of a discrete solution 



In this subsection, we derive, for Nc large enough, the existence of a unique 



local minimum of the discretized problem (77) in the neighborhood of 7rj^$°. 



Let 

Bn^ = {W^^ G n [7r^^*°]^ I < MwM^^vV^, < 1} , 
and f be the energy functional defined on Bn^ by 

SnAW") = {7T^^^° + 5(W^^-)7r^^$° + (103) 
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According to the fourth assertion of Lemma 4.2, the application 

defines a global map of Vj^ n M""^-^" such that C(0) = tt^^o. Therefore the 
minimizers of 

inf g v^AA (^04^ 

are in one-to-one correspondence with those of the minimization problem 

inf {SnAW^), e Bjv,} . (105) 

In a first stage, we prove that for Nc large enough, (105) has a unique solution 
in some neighborhood of 0. As a consequence (104) has a unique solution in 
the vicinity of 7r^$° (for Nc large enough). In a second stage, we make use of 
the unitary invariancc (69) to prove that for A^c large enough, (77) has a unique 
solution in the vicinity of 

Lemma 4.7 There exists r > and such that for all Nc > Nc , the func- 
tional has a unique critical point Wq" in the hall 

{w^^ G < n I WW^^Ui < r} . 

Besides, W^" is a local minimizer of (105) and we have the estimate 

W^^\\Hl< -^\\^^y -nnl- (106) 



Proof We infer from Lemma 4.6 that 
£nAW') = + - + «5(I^'^°)7r^,$° + 

+2a$o ((7r5^^$" - $°) + 5(Vr^=)7r^^$° + VT^^ {-k^^ - $°) + 5(iy^=)7r^^$° + W^^-) 
+R ((tt^^^*" - $°) + S{W'^^)TT^y + VK^'O 
= £;KS (^0) ^ 2a$o(I^^% VK~=) +4a$o (M/^% - $°)) 

+2a$o(7r^^^$° - $0,77^^$° - $0) +7ejv,(l^^=) 

where 

T^nAWn:) = 2a*o(5(VF^=)<$°,5(W^^0<<I>") 

+4a$o(5(W'^=)7r5^^$°, _ $0) + w^^^) 

+i?((7r^^^$° - $°) +5(iy^«)7r^^$° + ir^=) . 
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Thus, 

+7^JV,(W'^»)-7^JV,(0). (107) 
It follows from Lemma 4.6, (88) and the continuity of a$o on (H^iT))-^ that 



\ # # # 

for a constant C-ji > independent of Nc. Let us introduce for Nc > and 
r > the ball 

Bjv,(r) = e Vj^ n I a$o(ly^^ VF^=) < r2a*o(7r5^^$o - $°,7r5^^$° - $°)} 

We deduce from Lemma 4.4, that for all r > and all > N*, we have 



Let rg = 2(2C$o/c$j,)'''/^. For all r > ro, there exists A^c,r > N* such that 

VA^c>^cr, 9B«(r)cB« and VTy^<= G fr), ||T4^^'^||ffi < 1. 
Therefore, for all r > Tq and all N^. > A^c,r we have dBj^^{r) C Sjv^ and 

SnAW"^) > fiv.(0)+c*c||W^^=||^. -4C^o||W^^=||^.|K$°-$°||h_. 

y # # # 

+2|KcI>0 - cf>"r + ||,>1<j,0 _ cf>0|| ||^A'c||2 \ 



„2 



- ^^=(o) + ^^(^-^o)iK*°-*°ii^i 

As ||7r1^$°--<f>'^||^i goes to zero when Nc goes to infinity, we finally obtain that 
for all r > ro, there exists some ^ > A'^^ such that for all Nc > N'c^^, 

dBNAr)cBN, and VW^= e 9BjVe(r-), SnAW^") > ^nA^)- 
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This proves that for each A^c > -^c 2ro ' has a minimizer Wq" in the ball 
i3]Vc(2ro). In particular, 



\\W, 



Nc 

II H 



1 < 



(108) 



Let he a critical point of £n^ such that ^W-^ °||j;,2 < i. We denote by 



= W^' - W^^% 



As both Wq'' and W""" are critical points of £m^^ we have 
so that 

Using the expression (107) for £Nc^ we can rewrite this equality as 

a<,o(^T^^=,<5w^^«) = 6^^«^<^<5w^=) + c*o($^^,$]v^,<^<^^w^^=), 

where 

hlsiW^^W^^dW""^) = -a$o((5(V^f'0 -5(^^o'^=))7r^^<i>^(5'(^^i^=) •<5W^^=)7r^^^$o + 5W^^= 
-a$o(((5'(T^f=) - S'{W^')) ■ ^W^^=)7r^^$°, (tt^^^^ - $°) +5(W^o'^=)7r5^^$° + l^'^) 
-a$o((5'(Wi^-) • ^VF^=)7r5^^$°,(5VF^=) 

and 



i?'($k - $°) • ((.S'(Wo^=) • (5W^^=)7r^^$° + (5W^^=) 



-i?'($]v - • ((^'(VFi^^) • (5VF^=)7r^^$° + ^W^°) 



Using Lemma 4.3 and (108), we obtain that there exists (7$o (depending only 
on $°) and such that for all Nc> N^, 



< 



On the other hand, remarking that for all £ Ai and all S'^ E T^,Ai, 
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and introducing the path (^(i))te[o,i]5 drawn on the manifold A4 and connecting 
and , defined as 

^(t) = $0 + s{tW^- + (1 - t)Wo^=)7r5^^$° + iWf = + (1 - t)W^% 
we obtain 



c*o ($^^ , , ^^0^^% <5W^^= ) 



E^^ (*(0)) • *'(0) - E^^ (*(1)) • *'(1) 



a$o(*(0) - *'(0)) + a<E.o(*(l) - *'(!)) 
1 ri 1 



As = (V'i(i), • • • , ^M{t)) e for all t & [0, 1], we have for all 1 < i < Af 
and all t e [0,1], 



dx 



so that 



i=i -^r 

^ i=i -^r 

= -E^? / (V'.(i)-0°)C(i)-^i^^'"(*°)(*'(i),*'W). 

Consequently, 

Using Lemma 4.5, we obtain 



L 
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As 

= {S'{tW^^ + (1 - t)W^-) ■ (5VF^-)7r^$° + SW^-, 
= {S"{tW^' + (1 - t)W^'){dW^%SW^-))Tr^^^°, 

we obtain that there exists some constant C G M+ such that for N^. large enough, 
Thus, 

This proves that there exists a constant r > such that for all A^c large enough, 
WW^" Wh^^ < r implies 6W'^- = 0. Hence the result. □ 

As the mapping Bjv,(2ro) 5 W^- ^ n-^^^° + S{W^)w^^^° + W^- defines 
a local map of V^^ n 7W in the neighborhood of we obtain that 

$0^^ = ttM^o ^ 5(Wp^c)^Ai$o ^ ^rN, ^he unique local minimizer of 

inf {£KS($;vJ, $iv, e ViO^ n A^<*°} , 
in the vicinity of 7rj^$°. Besides, 

< C||n;v.*°-$°||Hi, 

for a constant C independent of A^c- We then have 

ll^n^,*o - IaaIIf < ||< - < c||n,v.$° - 

Let = [/jo^ _^o$^,, where C/jo^ = M|„ ^^(Mj^ _^oM|, ^J-V^. 

Then for each A^c > -^c2ro' ^at unique local minimizer of (77) in the 

set 



for some constant r*^ > independent of Ac, and it satisfies 

11*^. - < c||njv,$° - (109) 
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for some C e K+ independent of Nc- 

As ^% e M.^° , we can decompose ^% as 

$o,^ = $° + 50,^$o + W^^, (110) 
where = S{W^J and VT^^ G (note that M^]^^ ^ Vj^ in general). As 

\\s%J[f<\[w'nJ\Ii (111) 

and ||4>^ — $°||£fi goes to zero when Nc goes to infinity, we have, for Nc large 
enough, 

^IIW^^IIl^ < 11*^. - < 2||VF^JU|, (112) 

1\\W"^J\h^^ < 11*^. - < 2||W^1|U^. (113) 

The discrete solution satisfies the Euler equations 

where = p$o^ and where the AfxAf matrix A^^ is symmetric (but generally 
not diagonal). Of course, it follows from the invariancc property (69) that (77) 
has a local minimizer of the form U with U GU {Af) for which the Lagrange 
multiplier of the orthonormality constraints is a diagonal matrix. 

4.3 A priori error estimates 

We are now in position to derive a priori estimates for H^"^^ — ^"Hffj and 
(A^ — A*^), where we recall that Aq = diag(e5, • • • , e^). 

Using (2), (109) and the inverse inequality (49), we obtain for each s > 1 such 

/ X A/" 

that $° e fi?|t(r)) and each 1 < r < s, 

< CNr^<S>%^ - Hat, $"11^1 + IIHat,*" - $"11^. 

< cn:-' (||$^^ - <i>°||^. + - n^r^cooji^, ) + ||n^^$o - 

< CN-^'--^\\ILn^^°-^°\\hs^. (114) 

In particular, for s = 3 and r = 2, we obtain that ^% converges to $° in 
(H^iT))^, hence in (L^(r))^. 
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We then proceed as in (37) and remark that 

I / T/Coulomb /- _ „0\ 

+ X(^(^c + P°^.J-'^(Pc + /)j<^.</>,V. (115) 

As, from (110), 

we easily obtain, using the convergence of $^ to in (Jf^(r) D L^(r))-^, 

I|A°^.-A°I|f 0. (116) 

Nc-i-oo 



For W e (Z/^(r))^, we introduce the adjoint problem 

find "Ifw e ^""'^ such that 
VZe$0'^, a$(vI/vy,Z) = (iy,Z)i2 



(117) 



the solution of whom exists and is unique by the coercivity assumption (73). 
Clearly, 

\\^w\\H^^<C\\Wh.^. (118) 
In addition, it follows from standard elliptic regularity arguments that 

\\^w\\Hl<C\\Wh.^, 

yielding 

W'fw - UnMlI < CN-^W\\l2^ (119) 
W^w - Un^^wWh^^ < CN-'\\WU.^. (120) 
Denoting by \t' = \l/$o and using (110), we get 

N c 

= {^X. - 5^^, $")l| - 0^0(^1/, ^0,^$°) + a$o(vI/, $0,^ - $") 
= ($^^ - 5°,^$°)^2^ - a$o(*, 5°,^$°) + a$o(* - njv.M/, - $°) 
+a$o(njv,*,$^ -$°). (121) 
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Prom the definition (71), the last term in the above expression reads 

N N 

so that from the definition of the continuous and discrete eigenvalue problems 

N N . 

+4 E E(^ii,JVe - >^ij) / iVeniV, Vi- (122) 

i=i j=i -^r 
The definition of IIjv,, and the fact that * e yields 



which finally provides the estimate 

TV AT „ 

+ E E(^?i,iVe - A° ) / (<^° jv. - </'°)V'i. (123) 



Using Lemma 4.5, (109), (111) and (120), we infer 



+ ||A°,^-A°||p||$°,^-$°|U^j. (124) 

We thus obtain, using (116) and the above estimate, that asymptotically, when 
Nc goes to infinity, 

Reasoning as in (114), we obtain that for each s > 1 such that $° G [H^{T) j 
and each < r < s, there exists a constant C such that 

- < CN-(^-^^\\Um^^' - ci>0||«j. (125) 

To proceed further, we need to make an assumption on the regularity of the 
exchange-correlation potential. In the sequel, we assume that 

• either the function p ^ e^°^(p) is in CH([o,+oo)); 
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• or the function pc + is positive everywhere. As it is continuous on , 
this is equivalent to assuming that there exists a constant 77 > such that 
for all a; e M^, Pc{x) + p'^{x) > r]. 

It follows by standard elliptic regularity arguments that $° then is in (H^~^^^^~'^ (F))-^ 
for any e > 0, and we deduce from (125) that (78) holds true for all < s < 
m+ 1/2. 

Then, following the same lines as in the proof of (38), we obtain the estimates 



/ T/Coulomb / 0^ 



<C\\p%^-p'\\h-^, 



and 



J LDA de^^^ 
- (Pc + P°)-^(Pc+P°«))<^%>V 



valid for all < r < m — 3/2. Using these estimates in (115), we are lead to 
|A°-,^, - A0.| < C - ^%J\%.^ + \\p%^ - p°\\h-.) . 

Now, 



\\pk-p"\\H, 



sup 



Noticing that 



p%. -p' = J: K^J' - E = E«iv. - + ^i)' 

i=l i=l i=l 

we deduce 

\\p%^_-Ah-^<C\\K.-^"\\h-^^ (126) 
since converges, therefore is uniformly bounded in H^{r). Thus 

|A^^ -A°\<C (11$^^ - + C|K - $°||^-.) . (127) 

The derivation of estimates for \\^%^ — follows exactly the same lines 

as the derivation of the estimate: starting from the definition 

(l^,$o -$0)z.. 



\\^%.-^°\\H-r= sup 



# 



and remarking that the solution '^w to the adjoint problem (117) satisfies 
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we proceed as in (121) to get 



(W,$o,^-$0)i| = (W^,50,^$°)l| + WW^1)l| 



= {W, S%^^°)l2^ + a^oi^w, W°,J 

c # c c 




that yields 




The proof of (78) follows and then we get easily from (127) that 



- A°||f < CeTV-^^™-!-^). (130) 



Hence (79). Finally, (80) is a straightforward consequence of Lemma 4.6, (73), 
(96), and (109). 

4.4 Numerical results 

In order to evaluate the quality of the error bounds obtained in Theorem 4.1, 
we have performed numerical tests using the Abinit software [9] (freely available 
online, cf. http://www.abinit.org), whose main program allows one to find the 
total energy, charge density and electronic structure of systems (molecules and 
periodic solids) within Density Functional Theory (DFT), using pseudopoten- 
tials and a planewave basis. 

We have run simulation tests with the Hartree functional (i.e. with e^^^ = 
0), for which there is no numerical integration error. In this particular case, 
the problems (74) (solved by Abinit) and (76) (analyzed in Theorem 4.1) are 
identical. 

For TrouUier-Martins pseudopotentials, the parameter m in Theorem 4.1 is equal 
to 5. Therefore, we expect the following error bounds (as functions of the cut-off 



energy £e = i(^)') 




< 



< 



(131) 
(132) 



< 



i-4.5+e 



(133) 
(134) 



< 



1-4. 5+e 
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The first tests were performed with the Hydrogen molecule (H2). The nuclei 
were clamped at the points with cartesian coordinates ri = (—0.7; 0; 0) and r2 = 
(0.7; 0; 0) (in Bohrs). The simulation cell was a cube of side length L — 10 Bohrs. 
The so-obtained numerical errors are plotted in log-scales in Figures 1 and 2. 
The second series of tests were performed with the Nitrogen molecule (N2). 
The nuclei were clamped at positions ri — (—0.55; 0; 0) and r2 = (0.55; 0; 0) (in 
Angstroms), and the simulation cell was a cube of side length L = 6 Angstroms. 
The numerical errors for N2 are plotted in Figures 3, 4 and 5. The reference 
values for e° and /^^ for both H2 and N2 are those obtained for a cut-off 
energy equal to 500 Hartrees. 



Relative energy error curve as function of Ec (log/log scale) 
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Figure 1: Error on the energy as a function of for H2 
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Relative HI - Norm error curve as function of Ec (log/log scale) Relative L2 - Norm error curve as function of Ec (log/log scale) 




Figure 2: Errors on ||4>^^ -<I>0|l^i (left) and \\<^%-'P"\\l2 (right) and \e1 j,^-e°\ 
(bottom) as functions of Ec for H2 



These results are in good agreement with the a priori error estimates (131)-(134) 
for both the H2 and N2 molecules 
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Figure 3: Errors on ||$5v, ~ ^°\\h^ (left) and \\^%^ - (right) as functions 

of Ec for N2 " # " # 
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